from sympy import *
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import numpy as np

'''
这里的p是先假设其为一个常数，用数学的方法算出x的表达式，利用吸收半衰期进行p的求解
这里的u同样是先假设一个常数，用数学的方法求解出y的表达式，利用排除半衰期求解u
'''
p = 0.1386  # 药物向血液系统中的转移率和药量的比值
u = 0.1155  # 血液系统中的药物排除率和药量的比值
k = 1100  # 药剂总量


def dx_dt(x, t):
    return np.array(-p*x+0*t)


# 根据公式列出y的微分方程
def dy_dt(y, t):
    x = np.exp(-p*t)*k  # x关于t的函数，通过Sympy求出
    return np.array(p*x-u*y+0*t)

# 在适当时候进行急救后的微分方程
def dz_dt(z, t):
    p = 0.1386
    x = np.exp(-p * t) * 1100
    u = 0.2310
    return np.array(p * x - u * z)

# 求解x关于t的方程
t = Symbol('t')  # 时间
x = Function('x')(t)
end = dsolve(Eq(diff(x, t, 1), -0.1386*x))
print(end)


t = np.arange(0, 25, 1)  # 设置时间序列
# 利用odeint求解微分方程
a = odeint(dx_dt, 1100, t)
b = odeint(dy_dt, 0, t)
c = odeint(dz_dt, 0, t)
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
plt.plot(a,c='r')  # 对应x(t)
plt.plot(b,c='b')  # 对应y(t)
plt.plot(c, c='y')
plt.ylim(0,1200)
plt.grid(ls='--',color='k',alpha=0.5)
plt.show()
